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2.  Objectives 


The  following  were  the  principal  objectives  of  this  project. 

1.  Numerical  analysis  of  a  linearized  formulation  for  three-dimensional 
buckling  and  stress-stiffening  phenomena.  This  work  complemented  an 
implementation  on  the  commercial  program  STRESS  CHECK. 

2.  Investigation  of  hp  methods  for  computations  over  woven  composite 
materials. 

3.  Modeling  of  multistructures  via  hp  interface  mixed  methods  for  domain 
decomposition  and  concatenation.  This  work  was  originally  begun 
to  complement  an  implementation  in  the  commercial  program  MSC- 
NASTRAN. 

4.  Development  of  hp  methods  that  are  robust  and  locking-free  for  elas¬ 
ticity  problems. 


3.  Summary  of  Accomplishments/Status  of  Effort 


1.  The  investigation  of  robustness  for  hp  methods  for  the  linearized  buck¬ 
ling  formulation  implemented  in  STRESS  CHECK  was  completed.  It  was 
found  that  the  method  was  robust  except  in  two  main  cases:  (a)  when  the 
domain  is  not  ‘thin’  and  (b)  when  over-refinement  is  used  in  the  presence  of 
singularities.  It  was  established  that  otherwise,  one  can  be  assured  of  accu¬ 
rate  answers.  The  results  were  communicated  to  ESRD,  the  formulators  of 
STRESS  CHECK. 

2.  A  specialized  element  for  the  study  of  woven  composite  materials,  which 
incorporates  various  aspects  of  the  structure  of  the  material  was  formulated 
and  tested.  It  gave  excellent  results  in  predicting  the  macro-  and  micro¬ 
mechanical  response  of  the  material. 

3.  The  search  for  a  non-conforming  hp  method  that  can  be  used  for  do¬ 
main  decomposition  and  concatenation  was  completed  and  several  meth¬ 
ods  formulated  and  analyzed.  The  results  were  communicated  to  MacNeal- 
Schwendler,  the  formulators  of  MSC-NASTRAN. 

4.  A  locking-free  hp  mixed  formulation  that  is  robust  for  elasticity  problems 
even  in  the  presence  of  curved  boundaries  was  formulated  and  analyzed. 
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4.  Accomplishments 


Numerical  analysis  of  a  buckling  model 

The  area  of  failure  prediction,  in  the  context  of  e.g.  electronic  compo¬ 
nents,  ceramics,  laminated  composites,  etc,  has  been  of  long-standing  impor¬ 
tance  to  USAF  needs.  Experimental  determination  of  loads  that  can  cause 
buckling  and  other  failues  is  both  expensive  and  unreliable.  Mathematical 
models  used  in  current  engineering  practice  are,  on  the  other  hand,  based 
on  dimensionally  reduced  descriptions  of  an  elastic  body.  The  assumptions 
necessary  for  such  dimensional  reduction  to  hold  can  often  result  in  a  large, 
unknown  modeling  error  when  non-zero  initial  stress  states  are  present,  lead¬ 
ing  to  inaccurate  failure  predictions. 

A  linearized  model  for  buckling  and  stress-stiffening  has  been  imple¬ 
mented  in  the  hp  code  STRESSCHECK.  This  model  does  buckling  analysis 
for  the  fully  three-dimensional  problem  at  hand,  rather  than  some  asymp¬ 
totic  (dimensionally  reduced)  limit.  It  finds  the  smallest  positive  multiple 
A  of  an  existing  (pre-buckling)  stress  state  that  will  result  in  buckling.  The 
use  of  the  hp  method  enables  solutions  over  singular  domains  to  be  well  ap¬ 
proximated,  and  ensures  that  no  locking  takes  place  even  when  the  domain 
is  very  thin. 

However,  a  potentially  serious  danger  of  the  method  is  that  it  character¬ 
izes  A  as  the  lowest  positive  spectral  value  of  a  non-compact  operator  T.  Such 
non-compact  spectral  value  problems  can  be  notoriously  ill-behaved,  due  to 
the  presence  of  spurious  approximate  eigenvalues,  which  can  completely  pol¬ 
lute  the  results. 

The  results  we  have  obtained  in  this  project  establish  that  (1)  spurious 
eigenvalues  are  absent  for  problems  of  engineering  interest  for  which  the  pre¬ 
buckling  stress  (To  is  bounded  (2)  spurious  eigenvalues  are  always  present  for 
problems  where  (Tq  is  unbounded.  For  the  latter  case,  we  demonstrate  how 
the  reliability  of  the  computations  can  still  be  assessed,  using  the  eigenvectors 
(the  buckling  shapes).  We  describe  our  results  below. 

1.  Absence  of  spurious  eigenvalues  for  bounded  (Tq 

Suppose  Sb  denotes  the  Browder  spectrum  of  T,  i.e.  the  set  of  all  points 
in  the  spectrum  which  are  not  isolated  eigenvalues  of  finite  multiplicity.  Then 
polluting  eigenvalues  only  appear  in  the  interval  [minEB,maxEB]-  By  con¬ 
sidering  the  asymptotics  of  a  plate  of  thickness  d  ,  we  have  shown  that  the 
Browder  spectrum  lies  in  an  0{(f)  interval  around  the  origin,  so  that  spectral 
pollution  will  not  occur  for  the  eigenvalues  of  interest,  provided  the  domain 
is  thin.  (The  thin  case  is  the  case  of  engineering  interest.)  See  [5,7]. 

Again  considering  the  asymptotics  of  a  plate,  we  have  investigated  (with 
Monique  Dauge)  the  effect  of  various  initial  stress  states.  We  have  shown 
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that  the  only  case  where  spurious  eigenvalues  could  occur  is  the  one  where 
the  initial  stress  is  of  pure  bending  type.  This  case  is,  however,  uninteresting 
from  the  point  of  view  of  applications.  For  all  initial  stresses  that  occur 
in  applications  (these  can  be  decomposed  into  a  bending  and  a  non-zero 
membrane  stress),  spectral  pollution  does  not  occur  in  the  region  of  interest. 

Our  asymptotic  analysis  has  also  established  that  boundary  layers  occur¬ 
ring  in  the  solution  do  not  lead  to  adverse  affects  in  the  approximation  of 
the  desired  eigenvalues. 

Finally,  we  have  characterized  the  initial  stresses  that  lead  to  pollution- 
free  eigenvalue  approximations  when  parts  of  the  domain  are  not  thin. 

2.  Reliability  of  computations  for  unbounded  co 

For  domains  with  corners,  or  domains  with  a  change  in  the  type  of  bound¬ 
ary  conditions,  it  is  well-known  that  the  stresses  can  have  r~°‘  singularities. 
We  show  that  in  this  case,  the  Browder  Spectrum  will  no  longer  be  0{(JP),  but 
0{(P*M)  where  M  is  dependent  on  the  maximum  of  the  stresses.  Hence,  the 
Browder  Spectrum  may  be  large  enough  to  absorb  the  required  eigenvalue. 

In  applications,  one  is  interested  in  the  lowest  eigenvalue  of  a  family  re¬ 
lated  to  the  spectrum  of  a  limiting  problem  (the  d  =  0  case) .  We  characterize 
the  values  of  a  for  which  this  limiting  spectrum  is  still  well-defined,  and  well- 
behaved.  As  long  as  the  refinement  is  not  too  severe,  we  have  demonstrated 
that  we  can  still  recover  the  lowest  eigenvalue. 

In  Figure  1,  we  show  typical  buckling  shapes  associated  with  non-spurious 
modes  (Fig  1(a)  and  1(b))  and  spurious  modes  (Fig  1(c)),  for  a  plate  with 
a  5  degree  wedge  cut  out  of  it.  The  solution  has  a  strong  singularity  at  the 
wedge-tip.  It  is  seen  that  for  spurious  modes,  the  eigenfunction  is  almost 
constant  over  the  whole  domain,  except  for  the  region  surrounding  the  point 
of  singularity  (see  the  detail  of  an  ‘inner’  element  in  Fig  1(d)).  Figure  2  shows 
the  convergence  of  the  spurious  modes  when  an  hp  method  is  used  with  n 
layers  of  refinement  at  the  singularity.  We  see  that  the  spurious  eigenvalues 
converge  to  0  at  0(e~^"),  and  are  independent  of  d.  The  limiting  eigenval¬ 
ues  (the  ones  of  interest)  are  almost  independent  of  n  (they  have  already 
converged).  These  decrease  at  an  0{(f)  rate  as  d  is  decreased.  Hence,  they 
can  be  recovered  as  long  as  they  are  smaller  than  the  spurious  eigenvalues, 
as  shown  in  Figure  2. 

Two  important  facts  that  emerge  from  our  investigation  are:  (1)  The 
degree  of  localness  of  the  eigenfunctions  can  be  an  effective  tool  in  assessing 
the  reliability  of  the  computations  (2)  There  is  often  an  optimal  level  of  layer 
refinement,  over  which  spurious  eigenvalues  can  arise. 

The  above  algorithm  plays  an  essential  role  in  an  overall  hp  strategy 
for  numerical  computations.  Our  results  indicate  that  subject  to  the  above 
limitations  (which  are  not  serious  in  practice),  the  method  is  robust  for  en- 
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Figure  1:  Buckling  modes  for  thickness=0.1:  (a)  1st  mode,  global  (b)  2nd 
mode,  global  (c)  3rd  mode,  local  (spurious)  (d)  Detail  of  3rd  mode 

gineering  applications. 

Numerical  methods  for  woven  composites 

The  stated  research  objectives  of  AFOSR  place  significant  emphasis  on 
the  use  and  performance  prediction  of  highly  heterogeneous  composite  media, 
capable  of  sustained  use  in  adverse  loading  environments.  These  include,  in 
particular,  polymer  matrix  composites  (PMCs)  and  ceramic  matrix  compos¬ 
ites  (CMCs).  The  structure  of  such  materials  is  quite  complex,  with  bundles 
of  fibers  (the  tows)  woven  together  and  imbedded  in  a  fully  dense  surround¬ 
ing  material  (the  matrix),  with  large  voids  usually  present  in  CMCs.  This 
complexity  results  in  corresponding  problems  in  their  accurate  numerical 
resolution.  Previous  computational  efforts  employ  various  modeling  simplifi¬ 
cations,  with  consequent  repercussions  in  the  resolution  that  can  be  achieved. 

Our  goal  has  been  to  develop  an  hp  method  for  such  materials.  The 
starting  point  is  to  exploit  the  unit  cell  (Figure  3)  which  repeats  throughout 
the  material.  This  local  element  is  modeled  to  incorporate  the  geometric  and 
material  features  of  the  material  (with  care  given  to  accurately  modeling  the 
geometry  of  the  tows,  for  instance).  Previous  work  by  Prof.  Charalambides 
of  the  Mech.  Eng.  Dept  at  UMBC  has  involved  solving  a  series  of  canonical 
local  boundary  value  problems  on  the  local  cell,  to  predict  the  response  at  this 
level.  The  local  problems  are  first  reduced  to  2-d  models  by  a  Kirchhoff  type 
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Computed  eigenvatue 


Figure  2:  First,  second  and  spurious  e-values  for  various  thicknesses  d 


Figure  3:  The  plain  weave  composite  repeating  unit  cell 
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(a)  Micro-stxess  (Matrix) 


(c)  Micro-stress  (Matrix) 
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Figure  4:  The  micro-stress  Gxx  of  (a)  the  inter-tow  matrix  (b)  the  tows.  The 
micro-stress  Wyy  of  (c)  the  inter-tow  matrix  (d)  the  tows. 

hypothesis,  and  then  solved  by  a  p  (spectral)  method  to  avoid  subdividing 
this  cell.  It  has  been  empirically  shown  that  the  accuracy  is  comparable  to 
that  using  a  20-noded  brick  element. 

In  joint  work  with  Prof.  Charalambides  [6],  we  obtained  the  global  re¬ 
sponse  by  assembling  the  local  stiffness  matrices  and  solving  the  resulting 
problem.  Non-linear  behavior  is  incorporated  by  updating  the  material  prop¬ 
erties  selectively  over  the  elements,  based  on  the  computed  strains.  The  key 
factor  is  that  mesh  refinement  never  destroys  the  periodicity,  since  accuracy 
is  increased  only  by  increasing  the  polynomial  degree.  (The  initial  model  as¬ 
sumes  some  simplifications,  such  as  no  provision  for  delamination  or  slipping 
of  tows  being  included  at  this  stage.) 

We  have  formulated  and  implemented  this  specialized  element  and  ob¬ 
tained  computational  results  with  it.  In  particular,  we  have  compared  these 
results  to  those  obtained  when  calculations  are  performed  using  a  homog¬ 
enized  version  of  the  woven  composite  material.  The  macroscopic  results 
match  well  when  the  domain  has  a  large  number  of  cells  in  it.  The  ad¬ 
vantage  with  our  approach,  however,  is  that  in  addition  to  the  macroscopic 
response,  we  are  also  able  to  obtain  micro-strains  and  micro-stresses  in  vari¬ 
ous  individual  parts  of  the  material.  Such  results  are  not  possible  when  only 
homogenization  is  used. 
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Figure  4  demonstrates  this.  Here,  we  show  the  micro-stresses  in  both 
the  matrix  and  the  tow  part  of  the  composite.  The  problem  is  a  standard 
near-tip  crack  problem,  with  the  crack  tip  being  at  the  center  of  the  lower¬ 
most  boundary  of  the  square  shown.  The  material  is  a  CMC,  which  has 
a  central  circular  void.  It  is  clear  that  the  inter-tow  matrix  phase  and  the 
tows  phase  exhibit  different  stress  intensities,  consistent  with  their  intrinsic 
material  dissimilarities.  Moreover,  in  Figure  4(c)  (showing  the  micro-stress 
Wyy  in  the  matrix),  the  stress  concentration  in  the  vicinity  of  the  void  hole 
at  each  unit-cell  (represented  by  a  narrow  high-stress  band  perpendicular  to 
its  stress  direction)  is  clearly  captured.  Currently  available  analysis  tools 
that  we  are  aware  of  generally  use  equivalent  homogenized  materials,  and 
are  unable  to  capture  such  details.  (We  note  that  such  detail  is  needed  for 
failure  prediction  at  the  micro  level.) 

Our  future  goal  is  to  further  develop  this  element  by  validating  the  as¬ 
sumptions  that  have  been  used  in  its  design.  Also,  we  will  seek  ways  to 
improve  its  efficiency  (since  the  calculation  of  the  stiffness  matrices  currently 
requires  a  large  amount  of  effort).  We  are  also  exploring  methods  to  incor¬ 
porate  hp  techniques  of  capturing  crack-tip  behavior. 

Multistructure  modeling  via  hp  interface  methods 

Our  previous  work  has  shown  how  hp  discretization  of  a  complicated 
domain  can  be  distributed  over  several  independent  meshers  and  then  re¬ 
composed  by  mortaring  techniques.  In  [1,2],  we  formulate  new  mortar  meth¬ 
ods  that  have  optimal  stability  and  approximability  properties  with  respect 
to  hp  and  p-discretization.  Our  results  have  been  communicated  to  the  mak¬ 
ers  of  the  hp  code  MSC-NASTRAN,  where  a  similar  method  is  implemented. 
The  most-used  application  through  this  implementation  is  that  “library” 
meshes  on  previously  analyzed  components  are  now  being  directly  incorpo¬ 
rated  into  global  discretizations. 

Robust  hp  mixed  methods 

We  have  continued  our  investigation  of  hp  mixed  methods,  and  shown 
how  these  can  be  applied  to  a  problem  of  finite  elasticity  [3].  A  survey  of  hp 
methods  for  domains  like  plates,  shells,  etc,  has  been  prepared  for  a  special 
JCAM  issue  on  numerical  methods  [4]. 
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ement  for  the  study  of  woven  composites,”  Computational  Mechanics 
(in  press). 

7.  Manil  Suri  and  Monique  Dauge,  “Numerical  approximation  of  the  spec¬ 
tra  of  non-compact  operators  arising  in  buckling  problems.”  (preprint, 
2001). 


7,  Interactions/Transitions 

Presentations  at  meetings,  conferences  and  seminars. 

Mar  ’98  US  Naval  Academy,  Annapolis,  MD,  “The  hp  finite  element  method 
over  thin  domains”  (seminar) 
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Jun  ’98  “ICOSAHOM  ’98  (International  Conference  on  Spectral  and  High 
Order  Methods),  Herzeliya,  Israel,  June  22-26, 1998,  “Optimal  approx¬ 
imation  of  singularities  by  non-conforming  hp  finite  element  methods” 
(invited  minisymposium  talk) 

May  ’99  “SIAM  Annual  Meeting,”  Atlanta,  GA,  May  12-15,  1999,  “Non- 
conforming  hp  finite  elements:  Mortaring  techniques”  (invited  min¬ 
isymposium  talk) 

Jul,  ’99  “ICIAM  ’99,”  Edinburgh,  Scotland,  Jul  5-9,  1999,  “Optimal  con¬ 
vergence  rates  of  hp  mortar  finite  element  methods”  (invited  minisym¬ 
posium  talk) 

Jul,  ’99  “ICIAM  ’99,”  Edinburgh,  Scotland,  Jul  5-9,  1999,  “On  the  approxi¬ 
mation  of  spectra  for  buckling  problems”  (invited  minisymposium  talk) 

Oct,  ’99  “AMS  Sectional  Meeting,”  Austin,  TX,  Oct  8-10, 1999,  “hp  meth¬ 
ods  for  buckling  problems”  (invited  talk) 

Apr,  ’00  “Elastic  Shells:  Modeling,  Analysis  and  Numerics”  MSRI,  Berke¬ 
ley,  CA,  Apr  17-28,  2000,  “The  numerical  analysis  of  an  hp  algorithm 
for  the  approximation  of  buckling  problems”  (invited  plenary  talk) 

May,  ’00  “p  and  hp  Finite  Element  Methods:  Mathematics  and  Engineer¬ 
ing  Practice”  St.  Louis,  MO,  May  31-June  2,  2000,  Organizer  of 
Conference,  “Reliability  of  an  hp  algorithm  for  buckling  analysis” 

Jun,  ’00  “IASS  -  lACM  2000  Fourth  International  Colloquium  on  Compu¬ 
tation  of  Shell  and  Spatial  Structures”  Chania,  Crete,  June  4-7,  2000, 
“Reliability  of  an  hp  algorithm  for  buckling  analysis”  (invited  minisym¬ 
posium  talk) 

Apr,  ’01  “Conference  in  Honor  of  Jim  Greenberg’s  60th  Birthday”  Carnegie- 
Mellon  University,  Pittsburgh,  April  20-21,  2001,  “The  approximation 
of  the  spectra  of  non-compact  operators  arising  in  buckling  analysis” 
(plenary  talk) 

May,  ’01  Colloquium  talk.  Dept  of  Mathematics,  University  of  New  South 
Wales,  Sydney,  Australia,  May  16,  2001,  “The  approximation  of  the 
spectra  of  non-compact  operators  arising  in  buckling  analysis” 

May,  ’01  Colloquium  talk,  Tata  Institute  of  Fundamental  Research,  Mum¬ 
bai,  India,  May  31,  2001,  “Problems  in  finite  elements” 
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Jun,  ’01  “ICOSAHOM  ’01  (International  Conference  on  Spectral  and  High 
Order  Methods),  Uppsala,  Sweden,  June  11-15,  2001,  “Reliability  of 
an  hp  finite  element  method  for  buckling  analysis”  (invited  minisym¬ 
posium  talk) 

Transitions 

1.  Organized  conference  on  p  and  hp  methods  in  St.  Louis  with  support 
of  AFOSR.  Dominant  thrust  of  the  conference  was  the  dissemination  of 
AFOSR-supported  (and  other)  research  in  p/hp  finite  elements  by  mathe¬ 
maticians  and  engineers  to  the  industrial  community. 

2.  Research  on  non-linear  buckling  is  in  collaboration  with  Prof.  Barna 
Szabo  of  ESRD,  7750  Clayton  Road,  St.  Louis,  Missouri  63117  (phone:  314- 
645-1423).  Technology  transfer  is  occurring  through  the  implementation  of 
this  method  in  the  commercial  code  STRESS  CHECK. 

3.  As  reported  in  the  last  grant,  work  on  non-conforming  hp  methods  for 
multistructures  helped  provide  theoretical  validation  for  an  interface  method 
implemented  in  the  commercial  code  MSC-NASTRAN.  Contact:  Dr.  John 
Schiermeier,  The  MacNeal-Schwendler  Corporation,  815  Colorado  Blvd,  Los 
Angeles,  CA  90041-1777  (phone:  213-259-3832). 

8.  Inventions/patent  disclosures:  None. 

9.  Honors  and  Awards:  Appointed  to  editorial  board  of  SIAM  Journal 
on  Numerical  Analysis. 
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